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Abstract. We present a new formulation of the equations of motion used in smoothed 
particle hydrodynamics (SPH). The spatial resolution in SPH is determined by the smooth- 
ing length, h, and it has become common practice for each particle to be given its own 
adaptive smoothing length, hi. This has the advantage that the dynamic range that may 
be spatially resolved is greatly increased, but also has the drawback that additional terms 
which account for the variability of the smoothing lengths should be included in the par- 
ticle equations of motion in order to satisfy conservation requirements. We refer to these 
additional terms as terms. The difference between our approach and that of previ- 
ous implementations is that these V/i terms have now been included in the equations of 
motion, whereas they had previously been neglected. This is achieved by defining a func- 
tional form for the /ijS, that depends on inter-particle distances only, and then deriving 
the equations of motion using a Hamiltonian formalism. 

A number of test calculations, designed to compare the effects of including and ignoring 
these V/i terms, are presented. We find that the inclusion of the V/i terms has no detri- 
mental effects on the ability of SPH to model known problems, such as one-dimensional 
shock-tubes or the adiabatic collapse of cold gas spheres, with reasonable accuracy, and 
may in some cases lead to improvements in their qualitative outcome. For problems on 
which SPH has been shown to perform rather badly, because of poor energy conservation, 
we find that including the Vh terms results in a dramatic improvement. In particular, 
non-conservation of energy during a head on collision between identical polytropes can 
occur at the level of AE 10% when the V/i terms are neglected. When the Vh terms 
are included, however, then this error reduces to AE' ^ 0.8.% 

1 Introduction. 

Smoothed particle hydrodynamics was introduced in order to study problems in astro- 
physics involving the motion of compressible fluid masses of arbitrary geometry in three 
dimensions (Lucy 1977; Gingold &: Monaghan 1977). For a recent review see Monaghan 
(1992), and references therein. Particles are used to represent a sub-set of the fluid ele- 
ments that arise in the Lagrangian description of a fluid, and because spatial derivatives 
are calculated analytically from interpolation formulae, rather than on a grid, SPH is by 
its very nature adaptive and thus well suited to problems in which large density contrasts 
can occur, such as the fragmentation of self-gravitating gas clouds. 

The spatial resolution of SPH is determined by the smoothing length, h. In early 
formulations h was either taken to be constant, or else was a globally time-dependent 
function of the mean number density of particles in the system. More recently, however, 
it has become common for each particle to have its own time dependent smoothing length 
which adapts according to the local number density of particles. Full advantage is then 
taken of the Lagrangian nature of SPH, and the dynamic range in spatial resolution of the 
method is dramatically increased. The introduction of time varying smoothing lengths 
can, however, lead to serious problems with energy conservation in certain situations, as 
was highlighted in a recent paper by Hernquist (1993). The problem arises from the fact 
that the use of variable smoothing lengths means that additional terms should appear 
in the particle equations of motion. Up until recently, these additional terms have been 
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ignored, and as a consequence the particle equations of motion are no longer conservative. 

In a previous paper (Nelson & Papaloizou 1993), using a Hamiltonian formalism, we 
have derived a set of conservative equations for a barotropic fluid which include terms 
accounting for the variability of the smoothing lengths. In this paper these equations 
are generalised for the case of an adiabatic fluid in which the entropy on each particle is 
conserved. The treatment is then extended to include the effects of viscous dissipation. 

By means of simple numerical tests, a comparative study between the newly derived 
version of SPH and a more standard formulation of the method, including spatially and 
temporally varying smoothing lengths, is presented. An initial set of calculations is pre- 
sented, in order to compare the qualitative outcome of different versions of our code when 
it is used to model known problems, such as one-dimensional Riemann shock-tubes, and 
the adiabatic collapse of an initially isothermal gas sphere. The inclusion of additional 
terms in the particle equations of motion are not found to have any detrimental effects 
on the qualitative outcome of the calculations, and in some cases can lead to improve- 
ments as a result of more accurate energy conservation. A second set of test calculations 
is presented that is specifically designed to test the conservation properties of the SPH 
algorithm. It is found that errors in the conservation of energy or entropy can occur at the 
level of ~ 10%, when the additional V/i terms are neglected, during a head-on collision 
between identical polytropes. For the same calculation, it is found that these errors are 
reduced to 0.8% when the V/i terms are included. 

The basic hydrodynamical and thermodynamical equations of the problem are pre- 
sented in Section (2). The forms of these equations used in SPH methods arc derived 
and discussed in Section (3). In Section (4), calculations designed to compare the quali- 
tative performance of different versions of the SPH method are described. In Section (5) 
we describe and discuss the results from calculations designed to test the conservation 
properties of the different formulations. Calculations designed to test the effect of the V/i 
terms on the numerical diffusivity of SPH are presented in Section (6) . Finally, conclusions 
are drawn in Section (7), and we present a general discussion on issues arising from our 
calculations, and on the interpretation of the terms. 

2 Equations of motion. 

The equations of motion for a compressible fluid can be written 




(1) 



^ = _1VP-V$ + S, 
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denotes the convective derivative, p is the density, v the velocity, P the pressure, U the 
internal energy per unit mass, and 7i a term comprising all non-adiabatic heating and 
cooling rates per unit mass. This set of equations is supplemented by an equation of state, 
which can be written in the form 

P = P{U,p) (4) 

in which the pressure is written as a function of the internal energy per unit mass and the 
density, which may be adopted as the two thermodynamic variables defining the state of 
the fluid. 

Alternatively, we may introduce the entropy per unit mass, S, which satisfies the equation 

where T is the temperature. The quantities T, P, and U are then expressed as functions 
of the two thermodynamic variables p and S. We note the thermodynamic relation 

dS 

For the special case considered here of an ideal gas with constant specific heat capacity at 
constant volume Cy, and constant specific heat ratio 7, we have 

P=ij-l)pU, (6) 

P = Cyi^-l)pT, (7) 

and 

S = c,log(V-T). (8) 
In this case, it is also sometimes convenient to work in terms of the quantity 

P 

K = — , where 

/ s 

K = Cy{-f - l)exp — 
and is sometimes referred to as the entropic function. 



3 Numerical method. 

The above set of equations is solved numerically using SPH (Lucy 1977; Gingold &: Mon- 
aghan 1977). As described above, SPH uses particles to represent a sub-set of the fluid 
elements that arise in the Lagrangian description of a fluid. This set of particles is chosen 
in such a way as to ensure that the particles' mass distribution maps the fluid density, p. 
By calculating the trajectories of the particles according to the hydrodynamic equations, 
it is then possible to calculate the properties of the fluid at later times. 
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3.1 Hydrodynamical equations. 
3.1.1 Density. 

In order to reduce the statistical fluctuations arising from representing a fluid continuum 
by a finite set of particles, it is necessary to introduce a smoothing procedure. A smoothed 
quantity may be defined by the expression 

(/(r))= J f{v')W{\r-r'\,h)d^r' (9) 

where /(r) is an arbitrary function, 14/^ is a smoothing kernel, and h is the smoothing 
length. For a function /(r) defined only at a discrete set of points, a Monte Carlo repre- 
sentation of equation (|9|) may be written 

TV 

(/(r,)) = ^^/(r,m|r.-r,|,/i) (10) 
j=i PJ 

where nij is the mass of particle j, and pj is the density at the location of particle j. If 
the function /(r) is chosen to be the density, p(r), then we obtain 

N 

p{ri) = Y,m,W{\ri-rj\,h) (11) 
i=i 

where the angled brackets have been dropped for the sake of brevity, with the understand- 
ing that smoothed functions are used throughout in SPH. 

The spatial resolution of a calculation is determined by the smoothing length, h. Most 
implementations of the SPH algorithm employ temporally and spatially variable smoothing 
lengths in order to increase the dynamic range that can be spatially resolved, and hence 
take full advantage of the Lagrangian nature of the method. When variable smoothing 
lengths are employed each particle, i, will have an associated smoothing length hi. In the 
standard formulation of SPH, it is then necessary, in order to conserve momentum, to 
symmetrise expressions such as 

W{\ri-rj\,h) = Wij 

with respect to i and j. In Nelson & Papaloizou (1993) we set, following (Evrard 1988), 

h^h,, = ^^^. (12) 

Alternatively one can use the method suggested by Hernquist & Katz (1989) which sym- 
metrises using the replacement 

Wij ^ '^[W{\ri - Tjlhi) + W{\ri - rjlhj)]. (13) 

For comparison purposes, we have adopted the method of Hernquist & Katz (1989) in this 
paper. The smoothed density (|^) then becomes 

N . 

p{v,) = J2 ^jl^[Wi\r, -r,\, hi) + W{\v, - r, |, h,)] (14) 
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3.1.2 Equations of motion. 

We now have to form equations of motion for each particle. From equation (^) we obtain, 
denoting quantities associated with particle i by a subscript i, 

=Fp^i + FG,i + Smsc,i- (15) 

Here we have denoted the forces on particle i due to pressure, gravity and viscosity by 
Fp,i,FG,i, and S«sc,i respectively. 

3.1.3 A conservative system. 

In the early formulations of SPH a global smoothing length was used that could be either 
constant (Lucy 1977), or a function of time by relating it to the mean density of the 
system (Gingold & Monaghan 1977). In an attempt to increase the dynamic range that 
can be modelled, later formulations (Evrard 1988; Hernquist & Katz 1989) have replaced 
the global smoothing length by one with a spatial and temporal dependence, together with 
an appropriate symmetrisation with respect to particle pairs (see above). This has the 
consequence that, for barotropic or isentropic systems, the particle equations of motion 
are no longer conservative as they should be because the forces have not been derived 
from an appropriate potential function. It has been suggested in the literature (Gingold 
& Monaghan 1982; Evrard 1988) that effects arising from the variation of the /ijS in time 
and space should be small, although there is no firm evidence to show that these terms 
may be neglected. This concern is also true for a system in which a time varying, global 
smoothing length is employed, as shown recently for example by Hernquist (1993), but 
it is expected that, in general, the effects of non-conservation will be less severe. This 
is because the time change in the smoothing length, which drives the non-conservation 
of energy, is attenuated by the fact that it is dependent on the global properties of the 
system rather than on purely local ones. 

In his recent paper, Hernquist (1993) has shown that, under certain circumstances, errors 
at the ~ 10% level can occur in the conservation of either energy or entropy, depending 
on whether the entropy or energy equations are integrated, when modelling adiabatic 
flows. These errors were shown to arise as a direct consequence of the use of variable 
smoothing lengths. In view of this, it is desirable that a conservative formulation for the 
various terms in the particle equations of motion be derived that can incorporate variable 
smoothing lengths. This has already been given for the pressure force in a barotropic fluid 
in Nelson & Papaloizou (1993). Here we generalise the approach to the case of a general 
adiabatic flow in which there may not be a barotropic equation of state, but in which the 
entropy of each particle is conserved. We then extend the treatment to include viscosity. 

3.1.4 The pressure force. 

An inviscid fluid with a barotropic equation of state is a Hamiltonian system, so it is 
desirable, when modeling such a system, that the particle equations of motion in SPH 
form a Hamiltonian system also. The thermal contribution to the total energy can be 
written (see Nelson & Papaloizou, 1993) as 
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U = F(p)p(fr, where — — = 

Jy rin r.'i 



dp P 

This has the Monte Carlo representation 

N 

U = Y.mkF{pk). (16) 

k=l 

It is then natural to use C/ as a potential function for deriving the pressure force on particle 
i through 

dU 

The same formulation goes through for an adiabatic system for which the entropy of each 
particle is conserved. One simply replaces F{p) in the above by the internal energy per 
unit mass IA{S^ p) written as a function of the entropy and the density. The Monte Carlo 
representation for the pressure potential function is then 

AT 

U = ^mkUk{Sk,Pk), (17) 

k=l 

where Sk denotes the entropy per unit mass for particle k. We remark that the property 
of U that is needed is 

'dU\ P 



dp) s P^' 

For an adiabatic flow Sk is constant, so that in both the adiabatic and barotropic cases 



the pressure force on particle i is found, after differentiating either ( [iq ) or (17) and using 
(|l|), to be 



Pk\ O \ 



Fp,i = -EE"^'fc"^i \^\^.ilW^Vk-Y,Vhk)^W{\Yk-v,\,h,)\ (18) 



k=\j=\ 



Pk 



which may be written in the form 




+ 



dW{r,j,hj) 



hiconst 



dri 



dW{rkj,hk) dhk dW{rkj, hj) dhj 



dhk 



dri 



+ 



dhi 



dri 



(19) 



where we adopt the notation rkj = rk — rj. The first part of the above pressure force, which 
does not involve derivatives of the smoothing lengths, is of the usual symmetrised form 
(Hernquist & Katz 1989). The second part, which involves derivatives of the smoothing 
lengths, arises because of the spatial and temporal variability of the hiS. We shall refer to 
terms of this type as V/i terms. In order to evaluate the second part, we need to specify 
the way that the smoothing length depends on the particle coordinates. In order that 
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the system conserve linear and angular momentum as well as energy, we take the hk to 
be functions only of the absolute distances between particles. Here, we shall adopt the 
prescription given in Nelson &; Papaloizou (1993), namely 



/ N N 

k=Gk{Yl XknHkn{\rk - Tnl) 
\n=l / 



(20) 



where Q and H are arbitrary functions, and Xkn are arbitrary numbers, being zero when 
k = n. After differentiation of equation (20) with respect to rj, we obtain 



dhk 



N 



^kiQ'k X! ( XknH'k 



n=l 



''kn'i 



+ Q'kXkiH'k 



ki 1 



nk_ 

I rife I 



(21) 



where a prime denotes differentiation of a function with respect to its argument. By 



combining equations ([19D and (|2l]) we eventually obtain, for the pressure force on particle 
i, the expression 



P,i 



N 



dW{r,j,hi) 



+ 



dW{rij,hjj 



hiconst 
N 



1 V" Pi ^Pj] 9^i^ij^^i) n' rrl 



1 



N 



Pi 



^ ^k^s'kXkiH'ki^ 



N 



m 



k=l 

N N 



n=l 

dW{rkj,hk) 
dhk 



1 -A ^, rik ^ Pj dW{rkj,hk) 



k=l 



\rik\ 



J 1 



dhk 



(22) 



This form of the pressure force leads to a fully conservative set of equations because it has 
been derived from either of the potential functions ( p!^ ) or ( [TtD which are explicit functions 
of the inter-particle distances only. 



3.1.5 The addition of viscosity. 

If dissipative mechanisms are included, then we no longer have a system that is Hamilto- 
nian. However, it is possible to ensure that the net rate of energy input into the system 
is equal to the time rate of change of the Hamiltonian. The equations derived above for 
inviscid systems are thus modified so that the energy dissipated through viscosity goes 
into thermal energy in a consistent manner. For calculations in which artificial viscosity 
is included, we add an artificial viscous pressure term, Hjj, (Monaghan & Gingold 1983), 
in only the leading term of equation (22), not involving V/i terms, such that 



^ 2 „2 ^ „2 ^ *J • 



PI 



P 



8 



With this addition 

^P,i ^ j -|- Smsc,i' 

Where it has been used, we have adopted the artificial viscous pressure given by Monaghan 
& Gingold (1983) in the form 

Uij = ^ (-afiijCij + Pfifj) (23) 

Pij 

where the notation Aij = ^{Ai + Aj) has been used, c is the adiabatic sound speed, and 

7 . . if V • r - • < n 

P., = { '^rl + ri' iiv,,.r,,<U ^^^^ 

otherwise. 

Here, the notation Vjj = (vj — Wj) has been used, and rf' = prevents the denomi- 

nator from vanishing. 

3.1.6 The entropy equation. 

For a representation of the entropy equation (0), we write, for particle i. 



dSi 



§ = n.- (25) 

„. at 

HI 



This may be incorporated into a conservative scheme (see below). 
3.1.7 Conservation of energy. 

In the case of either a strictly adiabatic or barotropic inviscid flow, the total energy, E is 
conserved. This may be written in the form 

E = Tk + u + n. 

Here 

1 ^ 

Tk = ■7:^rni\-Vi\^ 
^ i=i 

is the kinetic energy, U is the potential function representing the thermal energy which 
we used above to derive the pressure force, and represents the gravitational energy. 
Provided that Q is expressed as a function only of the distances between particles, it may 
be used to derive the gravitational force. 



3.1.8 The effect of viscosity. 

When viscous terms are included in the equations of motion, and we allow the entropy Si 
associated with particle i to vary with time, we may derive an equation for E by taking 
the scalar product of Vj with the equation of motion ( |l5|) and summing over i, so obtaining 
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dE -A dUi dSi 1 \^ \^ TT ^ij ' ^ij 

—— = > mi-— — — — - > > rriimAVij — r — x 
dt ^ 'dS, dt ^tij^i l^^^l 

dW{vij,hi) 
d\T\ij 

Here, the second term, involving the double summation, represents the viscous dissipation, 
which is necessarily positive definite. 
If we write equation in the form 



+ 



dW{Yij,hj, 



const 



d\r\. 



(26) 



dUj dSj 
dSi dt 



N 



4 ^ 



dW{rij,hi 



d\r\. 



+ 



dWirij,h,) 



hiconst 



d\r\ 



(27) 

where the first term on the right hand side represents the heating rate per unit mass due 
to viscosity and Ci represents the net non- viscous heating rate per unit mass, then we see 
that equation ( p^ ) becomes 

dE ^ ^ 

Thus the rate of increase of the total energy is equal to the net non viscous heating rate 
as required. 

The evolution of the entropic function, K, defined in Section (2), is then given by the 
expression 



dKjSi) ^ 
dt 

Using the second law of thermodynamics 



(7 - m^p. 
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(28) 



dUi dUi dSi Pi dpi 
~dt " 'dSl~dt ^ '^i~dt' 



we may derive an equation for the rate of change of the internal energy per unit mass for 
particle i. Using ( p7|) and differentiating equation this may be written in the form 



dUi 
~dt 



ni + 



1 r\ 

2 



N 



dW{r,,,h,) 



9|r 



+ 



1 n 



N 



N 



N 



dhi 

dW{vij,hj) I 



n=l 
N 



+ 

hiConst 

^in ' Tin 

I Tin I 



5|r 



dhj 



yj 2^ Xjn^jn , . I 



(29) 



n=l 



Equation (29) may thus be used as an equivalent alternative to the entropy equation (25). 
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At this point, we note that the standard method of deriving the SPH equations of motion 
- i.e. by multiplying the continuum equations by W, integrating over the solution domain, 
and integrating the resulting expressions by parts (see Hernquist & Katz 1989) - will not 
result in conservative equations of motion when spatially variable smoothing lengths are 
employed, even if the V/i terms that arise from such an approach are included. The reason 
for this is that, in the standard method, the pressure forces experienced by the particles 
can be expressed in the form of a sum over pairwise interactions. However, when the 
smoothing lengths are expressed in terms of inter-particle distances, it is apparent that 
the pressure force on a particle i can no longer be expressed in the form of a sum of simple 
pairwise interactions, but simultaneously involves communication with other particles in 
the system - namely all the particles which contribute to hi , and also those particles in the 
system that receive a contribution from particle i to the calculation of their own smoothing 
lengths. 

In principle, it is not necessary to employ either of the symmetrisation procedures de- 
scribed by equations (^) or (jl^) when the equations of motion are derived from a Hamil- 
tonian. When the Hamiltonian is formed using either of the so-called gather or scatter 
formulations (for definition of these terms see Hernquist & Katz 1989), the resulting equa- 
tions of motion are symmetric and will conserve both momentum and energy. However, 
the form in which they appear does not allow the artificial viscosity in equation (23) to 
be incorporated in the usual manner. 

3.1.9 The smoothing kernel. 

The smoothing kernel used is that proposed by Monaghan & Lattanzio (1985), and takes 
the form 



W{u,h) 




0<u/h<l 

l<u/h<2 
u/h > 2 



(30) 



where n = |rj — rj\. 

3.1.10 Smoothing lengths. 

There are two basic requirements that must be met by the method chosen to calculate 
the individual particle smoothing lengths. Firstly, the smoothing lengths must be explicit 
functions of the distances between particles, and secondly, the number of particles with 
which a given particle interacts should be roughly constant. We use a method similar to 
that suggested by Hernquist & Katz (1989), in order to keep the number of neighbours 
within 2hi,J\fTOL, constant. Having obtained this list of nearest neighbours, hi is then 
defined to be 

h^ = -\ri - r,,nax\ (31) 

where rimax is the position vector of particle i's most distant nearest neighbour. In terms 
of our general functional form for hi given by equation (pO|), the functions Q{x) = x 
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Xin — 'Jn,imax 

becomes 



, and H{x) = \x. From equation (p2|), the pressure force on particle i then 



1 ^ 

-Y 



rrnmj -9 + -5- 



3r 



+ 



dW{rij,hj] 



hiconst 



di 



TV 



Pi , 9VF(rij-,/i,) ri 

P: 



EPk tr ^kmax k 

^k 2" * kmax' ' / 



fc=l 

TV 



l^kmax k 



dhi _ , 

^ dW{rkj,hk) 



1 r ^kmax k 

T / ^ ^k'J'i kmax' 



k=l 



^kmax k 



N 

E 



1 J2 



dhk 



(32) 



and similarly for the equation governing the evolution of the specific thermal energy we 
obtain from (p9|) the expression 



dUj 
~dt 



dW{Yij, hi) Vi 




9|r 



+ 



hi const 



d\r. 



N 

1 jmax 
^j— 



jmax i^ij 7 f^j J 



J jmax\ 



(33) 



We note that, in order for the potential function given by equation (^) to be a continuous 
function of its arguments, the functional form adopted for the hiS must be continuous. 
Consequently, the number of nearest neighbours, Mtol, that determines the value of each 
hi should be strictly constant throughout the calculation, and particle imax in equation 
(|3l|) must always be the 'A/roLth' nearest neighbour. This point is discussed further in 
Section (5). 

We also note that the functional form (^) used to calculate the smoothing lengths 
introduces additional statistical fluctuations into the pressure forces. From the momentum 
equation ( |3^ ) it is apparent that when particle z's nearest neighbour, imax, changes its 
identity as the system evolves in time, fluctuations will arise in the pressure forces ex- 
perienced by particle i. Technically speaking, the potential function defined by equation 
(|T^ is itself continuous, but it has a discontinuity in its first derivative, the magnitude of 
which is determined by the V/i terms. An alternative method for calculating the smooth- 
ing lengths, which reduces these fluctuations, is to first find the Ntol nearest neighbours 
and to then calculate hi from this list of nearest neighbours according to 



1 



N 



N 



f"-"- n=l 



far ^ 



(34) 
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where the summation is over particle i's Nfar most distant nearest neighbours. In par- 
ticular, we have experimented with using Nfar = 6 so that the value of 2/ij then be- 
comes equal to the mean distance to particle Vs six most distant nearest neighbours. 
In terms of the general functional form for hi given by equation (^), the functions 
G{x) = x^x = ifar, H{x) = x/2Nfar, where each ifar represents one of particle i's 
Nfar most distant nearest neighbours. 



3.1.11 Gravitational forces. 

An implementation of the Barnes-Hut hierarchical tree algorithm (Barnes & Hut 1986; 
Hernquist 1987) was used in the calculations for which the computation of the force Fq,!, 
due to self-gravity was required. The interparticle potential was softened by the method 
of spline softening (Gingold & Monaghan 1977), with a constant value of the softening 
length used throughout. When calculating the gravitational force on a given particle i, 
an opening angle of = 0.6 (where 9 = s/d; s =size of cell, (i=distance to centre of mass 
of cell), was used for cells whose centre of mass lay outside of the softening radius, with 
forces being calculated up to quadrupole order. For cells lying within the softening radius, 
an opening angle of ^ = 0.3 was employed, and forces were computed up to monopole 
order only. 

This method of calculation gives an approximation to the gravitational potential in 
such a way that it is no longer strictly a sum of contributions dependent only on the 
distances between particles, and this will lead to some small (AE < 1%) error in the total 
energy conservation. 



3.1.12 Time integration. 

The standard second order leap-frog scheme, with the modification proposed by Hernquist 
Sz Katz (1989) for estimating time-centred velocities in the viscous pressure term (p3|), 
the energy equation (p9[), and the entropy equation (27), was used. The time step was 
determined by the expression 

6t = QMiui ^ ■ (35) 

Ci + l.2[aci + pmaxjlfiijl) 

where the factor Q is a numerical constant usually taken to be in the range Q ~ 0.1 — 0.3. 



4 Numerical tests. 

A number of test calculations have been performed using different versions of our SPH code 
in order to compare their ability to reproduce the solutions for problems that have either 
been studied using independent numerical techniques, or which have known analytical 
solutions. The main purpose of this paper, however, is to study the effect of including the 
V/i terms on the conservation of energy, and the emphasis will be on those calculations 
designed to test conservation properties. Tests have also been performed to explore the 
possibility that the inclusion of V/i terms in the manner described above may lead to 
increased numerical diffusivity. 
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4.1 Shock-tube calculations. 



Extensive one-dimensional Riemann shock-tube calculations (see Sod 1978; Monaghan &; 
Gingold 1983), have been performed. The initial conditions were: 

p = 1, = 1, f = 0, for X < 0; , . 

p = 0.25, P = 0.1795, = 0, for X > 0. ^ ' 

Here x is a cartesian coordinate and the calculations were for an ideal gas with specific heat 



ratio 7 = 1.4. The equations integrated were the equation of motion (15) using (^) and 
the entropy equation (p^) with Ci = 0. Plots of the density, pressure, velocity field, and 
entropic function K, are shown in Figs, (l.a - l.d) for the case when the V/i terms were 
neglected, and in Figs. (2. a - 2.d) for the case when the V/i terms were included. In both 
cases shown, 400 equal mass particles were distributed in the range —0.6 < x < 0.6 so as to 



satisfy equation (36). The smoothing lengths were allowed to vary in both space and time, 
with the number of nearest neighbours being fixed at Mtol = 4, using a renormalised, one- 
dimensional version of the smoothing kernel (|3^). This corresponds to ~ 32 neighbours for 
a spherical kernel in three-dimensions. For the calculations presented here, the artificial 
viscosity parameters in equation (^) were chosen to be a = l.,/3 = l.^rf = 0.01, and the 
time-step was calculated using equation ( p5[ ) with Q = 0.2. 

The results shown in both Figs. (1.) and (2.) are almost identical in form, and show 
a satisfactory agreement with the analytical solution. The shock front located between 
X = 0.2 — 0.25 is broadened over a range in x of « 2/ij, with the contact discontinuity, 
positioned at x ~ 0.1 having a similar width. The slight overshoot in the density and 
pressure at x « —0.2 appears to be due to the fact that the values of the smoothing 
lengths become smaller at this point in order to satisfy the constraint that the number 
of nearest neighbours be constant. This effect was seen to decrease when the number of 
nearest neighbours Mtol was increased, though this also had the effect of broadening the 
shock front. 

4.2 Gravitational collapse of cold gas sphere. 

The adiabatic collapse of an initially isothermal gas sphere has been calculated using 
different versions of our code in order to test the qualitative effects of including the V/i 
terms in a three-dimensional calculation. The initial conditions were chosen to match those 
of Evrard (1988) and Hernquist &; Katz (1989) so as to facilitate comparison with their 
P3MSPH and TREESPH codes respectively, and with a finite-difference code (Thomas 
1987). 

Following Evrard (1988), the initial state is an ideal gas sphere of radius R and mass Mt-, 
with a density profile given by 

where r is the usual spherical polar radius. The initial conditions were set up by placing 
the particles on a uniform grid which was then stretched to achieve the required density 
distribution. Initially, the gas sphere was isothermal, with specific thermal energy u = 
0.05GMt / R- The ratio of the specific heats is 7 = 5/3. 
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As in Evrard (1988), the scales used in presenting the results are: density, = ?>Mt /^ttB?] 
energy, = GMt/R; pressure, P=k = /0*u*; and velocity, v^, = {GMt / RY^"^ ■ Time is shown 
in units of the free-fall time at the initial outer radius, t^, = {tt"^ /SG)^^"^ R^^'^Mj,^^'^ . The 
equations integrated were the equation of motion (15) using (|3^ , and the entropy equation 
(27) with Ci = 0. The results for the calculation in which the V/i terms were neglected 
are shown in Figs. (3. - 6.), and those for which the V/i terms were included are shown 
in Figs. (7. - 10.). The gas sphere begins to collapse due to the low value of the internal 
energy, and as it does so the ratio of the thermal to the gravitational energy increases 
since 7 > 4/3. At later times, a central bounce occurs and a shock wave propagates 
outwards through the collapsing sphere, with the system converting most of its kinetic 
energy into thermal energy between times t ~ O.Stj, and t ~ 1.2t*. Eventually, the sphere 
tends towards an approximate equilibrium configuration with U ~ — J7/2. 

Both of the calculations displayed in Figs. (3. - 6.) and Figs. (7. - 10.) show excellent 
agreement with the SPH calculations of Evrard (1988), and those of Hernquist & Katz 
(1989). Particularly encouraging is the good agreement found between our calculations 
and the results obtained by Evrard (1988) using a one-dimensional finite difference code 
containing 250 zones (Thomas 1987). The strength and position of the shock at time 
t = 0.8t=K is well represented by our calculations, and the thermal energy profile behind the 
shock shows the slow rise with radius displayed by the one-dimensional model. Similarly, 
the velocity profiles from our calculations, shown in Fig. (6.) and Fig. (10.), indicate 
good general agreement with those of the one-dimensional calculation, though there is a 
larger scatter in these results than in those for the density and pressure caused largely by 
non-synchronous evolution of the system at different angular positions around the sphere. 

When comparing the results of our two SPH calculations, small differences may be 
observed. Firstly, the results for the case in which the V/i terms were included appear 
to suffer from a slightly larger scatter than the results from the calculation in which 
the V/i terms were neglected. This 'noise' appears to be due to the functional form, 
(|3l|), adopted for calculating the his, which, as described in Section (3.1.10), introduces 
increased statistical fluctuations into the pressure forces. However, the overall effect of 
these increased fluctuations appears to be small. Closer inspection of the results shows 
that including the V/i terms actually leads to improvements in the performance of the 
code due to energy being conserved more accurately. For example, at time t = 1.2t^, 
the sharp peak and drop in the thermal energy located at r « 0.6R in Evrard's one- 
dimensional calculation is more apparent in our results if V/i terms are included than if 
they are neglected. Consequently, the corresponding feature in the pressure profile is also 
more accurately represented. Energy was conserved to a level of ~ 1.0% when V/i terms 
were included, and to a level of « 4.3% when neglected. 

The calculation was repeated using the energy equation (33) with the V/i terms being 
both included and neglected, respectively. As expected, the qualitative outcome of these 
calculations were almost identical to their equivalent calculations displayed in Figs. (3. - 
10.) in which the entropy equation ( p7| ) was integrated. If the V/i terms are neglected, 
then no apparent qualitative improvement in the results is obtained by using the thermal 
energy equation (^) - which leads to conservation of total energy but not entropy, instead 
of the entropy equation (^) - which leads to the conservation of entropy but not energy. 
The calculation was also performed using the functional form for the /ijS given by equation 
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(|3^ ) with Nfar = 6, and Mtol taking the value 50 ± 1. In this case the V/i terms were 
included, and the entropy equation (|^) was integrated. The results were almost identical 
to those shown in Figs. (7.) - (10.), except that the level of scatter was reduced slightly. 
Energy was conserved to a level of ~ 0.6%. 

Obviously, the inclusion of V/i terms does not alter the ability of SPH to reproduce the 
salient features present in the problems studied thus far, and can in fact improve the 
qualitative outcome of the calculations. 



5 Conservation of energy. 

Two different sets of test calculations were performed in order to test the energy con- 
servation properties of different versions of our code. These were collision calculations 
between identical polytropes, including the effects of self-gravity and artificial viscosity, 
and the free expansion of non self-gravitating polytropes in which the effects of viscosity 
were ignored. 



5.1 Colliding polytropes. 

Numerical calculations of collisions between identical n = 3/2 polytropes have been per- 
formed. Equilibrium polytropes, each consisting of = 2085 particles, were set-up by 
evolving an initially uniform gaseous sphere in the presence of damping until a stationary 
state was achieved (Lucy 1977; Goodman & Hernquist 1991). Rather than being allowed 
to collide after free falling from a distance, the two identical polytropes were placed adja- 
cent to one another with their surfaces just in contact. The polytropes were then given the 
free-fall velocities that they would have if they had started from a state of rest, with an 
initial separation between their centres equal to six times their diameter. This procedure 
was followed in order to reduce the computational costs. 

Collision calculations for six different versions of our code were performed, the details of 
which are summarised in Table 1. In these calculations the equations integrated were the 



equation of motion ( [iq ) using (|32|), and either the entropy equation (27) or the internal 
energy equation (^3|). In all cases Ci = 0. A typical time sequence is shown in Fig. (11.), 
and the evolution of the total energy for a number of the calculations described in Table 
1. are shown in Fig. (12.). For all the calculations presented, a value of Q = 0.2 was 
used in equation ( ^5| ) to determine the time step, and the artificial viscosity parameters 
in equation ( [2^ ) were a = 0.5, /3 = 1. 

As the polytropes collide, a shock forms at the interface converting most of the kinetic 
energy into thermal energy. A period of rapid re-expansion then occurs due to the presence 
of large pressure gradients in the fluid, but because the system remains bound, the resulting 
single polytrope undergoes quadrupole oscillations which damp out as the system tends 
towards an equilibrium state. 



5.1.1 V/i terms neglected. 

From the results for the evolution of the total energy shown in Fig. (12.), it is obvious that 
large differences exist between the results for the two calculations in which the V/i terms 
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were neglected. It appears that the total energy is conserved with reasonable accuracy 
(at the ~ 2% level) if the energy equation (^) is integrated, whereas errors at the level 
of ~ 10% are incurred if the entropy equation ( p7| ) is integrated. These errors largely 
occur between the times t = 0.1 and t = 0.3 during the initial compression and rapid 
re-expansion phase. However, the apparent improvement in the conservation of energy 
obtained by integrating the energy equation is somewhat illusory since entropy is no longer 
conserved (Hernquist 1993). This is because neglecting the V/i terms in both ( |32|) and 
( p^ is consistent with conservation of total energy but is inconsistent with the correct 
equation for the entropy. We have performed further calculations that show that the 
errors incurred in the conservation of the total energy when integrating equation ( |33[ ) are 
largely due to the use of velocities that are not properly time-centred and from the use 
of a tree-code to compute the self-gravity of the fluid. Also, calculation C4 in Table 1. 
was repeated using a value for Q in equation ( |35|) of Q = 0.03. Violations in energy 
conservation still occurred at the level of 9.75%, implying that non-conservation is not 
induced by time-stepping errors, but is instead due to a physical change of the energy in 
the system introduced by the time variation of the smoothing lengths. 



5.1.2 V/i terms included. 

It is apparent from Fig. (12.) that the inclusion of V/i terms results in a dramatic 



improvement in the conservation of total energy. The entropy equation (27) was integrated 
for all runs displayed in Fig. (12.) in which the V/i terms were included. As mentioned 
in Section (3.1.10), by allowing the value of Mtol to vary from a single, fixed value, 
temporal discontinuities are introduced into the potential energy of the system, and energy 
is no longer exactly conserved. This point is illustrated by our results, which show the 
conservation of total energy slowly degrade as the value of Mto l is allowed to vary from 
J^TOL = 50 to Mtol = 50 ±1, and Mtol = 50±2 respectively. However, even with Mtol 
varying by ±2, the errors incurred are still relatively small (~ 3%), and for less violent 
scenarios are almost negligible. Computational advantage is gained by allowing Mtol 
to vary slightly about some fixed value rather than maintaining constancy. Any method 
used to estimate the /ijS at the beginning of a time step has a greater chance of satisfying 
the demand that the number of nearest neighbours lies within a small range about some 
value, rather than being equal to a single, specific value. Therefore, the estimated hiS 
have to be adjusted less often and the execution speed of the code is increased. Errors 
incurred in the conservation of energy for the run in which Mtol took a single value 
(i.e. Mtol = 50), were due to residual discretisation errors, and the use of a tree-code to 
compute gravitational forces. 

The collision calculation was also performed by including the V/i terms and integrating 
the thermal energy equation (33). Although we do not display the results from this 



calculation in Fig. (12.), because integrating either of the equations (|27D or ( |33|) is formally 
equivalent when V/i terms are included, we found that the evolution of the total energy 
was very similar to the case in which equation ( p^ ) was integrated and the Vh terms were 
neglected, as should be expected. The deviation from perfect energy conservation in this 
case is again due to the use of velocities that were not properly time-centred in equation 
and the use of a tree-code. 
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5.1.3 Timing tests. 



The CPU time per time step was monitored for a number of collision calculations in order 
to estimate the increase in computation time arising from the inclusion of the V/i terms. 
The results are presented in Table 2. The two parameters that were varied were the 
number of nearest neighbours within 2/ij, Mtol, and the functional form used to calculate 
the smoothing lengths, which was either equation ( [31] ) or equation (Q). A total of 4170 
particles were used in each of the calculations. As mentioned in Section (5.1.2), allowing 
the value of Mtol to vary about some value, 50 say, leads to a more efficient code, since 
less time is spent trying to iterate 2hi towards a value such that it exactly encloses the 50*^ 
nearest neighbour. It is apparent from Table 3 that the inclusion of the extra V/i terms 
leads to only a small ( < 10%) increase in the computational overhead when equation (31) 
is used to calculate the smoothing lengths. When equation (|^) is used, with Njar = 6, 
then the computational expense is greater, due to it being necessary to find the positions of 
the six most distant nearest neighbours from which the hiS are then calculated. However, 
a rather naive sorting algorithm was used to perform this task, so we expect to be able 
to decrease the computational requirements with a more sophisticated approach in future 
calculations. All calculations shown in Table 3 were performed on a SUN IPX workstation, 
operating at 33MHz with 16 Mbytes of core storage. 



5.2 Expanding polytropes. 

In order to test the different forms of the SPH algorithm under less extreme conditions, 
and in the absence of gravitational and viscous effects, a series of test calculations were 
performed on the free expansion of gaseous polytropes. 

An n = 3/2 equilibrium polytrope, consisting of = 4945 particles, was set up in the 
manner previously described in Section (5.1). At time t = 0, the self-gravity was switched 
off. As expected, the absence of pressure at the surface boundary causes a rarefaction 
wave to propagate through the polytrope which expands freely into space, with the ther- 
mal energy being converted into the kinetic energy of the fluid's bulk motion. Artificial 
viscosity was neglected, though relation (|3^ ) was used to calculate the time step, with 
a = 0.5, /? = 1., Q = 0.05. 

A number of free-expansion calculations were performed, with the run parameters being 
summarised in Table 3. The variation of the total energy in the system as a function of 
time is plotted in Fig. (13.) for a number of the runs described in Table 3. 
Although the free expansion of a non self-gravitating polytrope is a considerably less 
violent phenomenon than the collision between two identical polytropes, a similar trend is 
seen to occur in the results for the conservation of total energy, though the discrepancies 
are less marked. 



5.2.1 V/i terms neglected. 

From Fig. (13.) it is obvious that once again the use of the energy equation ( |33|) rather 
than the entropy equation (^) results in energy being conserved to a higher degree of 
accuracy when V/i terms are not included. If the energy equation is integrated then errors 
incurred in the conservation of the total energy are small, and mainly arise from the 
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fact that the velocities used in equation ( pS] ) are not correctly time-centred. In this case, 
however, it can be shown that entropy is not conserved. Errors at the level of ~ 1.5% arise 
if the entropy equation (^) is used instead of p^). This is because the total energy is not 
conserved if time varying smoothing lengths are used and the V/i terms are neglected when 
modeling strictly adiabatic flows. Thus, errors of comparable magnitude are incurred in 
the conservation of different fundamental properties of the fluid when variable smoothing 
lengths are employed, depending on whether one chooses to integrate the entropy equation 



(27) or the energy equation p^). We note that the time period during which most of the 
energy was lost from the system in calculation E3 (see Table 3.), i.e. t = 0.0 — t = 0.2, 
corresponds to ~ 1200 time steps in our simulations, and represents the point when more 
than 99.5% of the initial thermal energy in the system has been converted in to kinetic 
energy - implying that non-conservation is not largely due to time stepping errors. 



5.2.2 V/i terms included. 

It is apparent from Fig. (13.) that the inclusion of the V/i terms results in much improved 
energy conservation. If the value of Ntol is kept constant, and the entropy equation is 
used, then energy is conserved almost exactly, with the remaining deviation being due to 
residual discretisation errors. 11 Mtol is allowed to vary slightly about a fixed value (i.e. 
J^TOL = 50 lb 2), then small errors (~ 0.15%) occur, which are substantially smaller than 
the errors resulting from integrating the entropy equation and neglecting the V/i terms. 
When equation ( |3^ ) is integrated instead of equation (p7|), then small errors, the evolution 
of which are similar to those from calculation in Table 3., are incurred (see run i?5 in 
Table 3.). We do not present the results from this calculation in Fig. (13.) for similar 
reasons to those given in Section (5.1.2). 



6 Statistical fluctuations and numerical diffusivity. 

This last series of test calculations is concerned with the increased statistical fluctuations 
resulting from the inclusion of the V/i terms, and in particular the functional form, (|3lD, 
adopted for calculating the /ijS. In order to examine the importance of this effect, two test 
calculations were carried out. 



6.1 Pressure force fluctuations. 

Firstly, the pressure forces on a random sub-set of the particles were monitored throughout 
the free expansion of an n = 3/2 polytrope. An illustrative selection from these results 
are presented, that represent the forces experienced by the particle as a function of time 
in the x, y, and z directions respectively. Fig. (14. a) shows the forces experienced during 
a calculation for which the V/i terms were neglected, with the entropy equation being 
integrated, and Fig. (14.b) shows the forces on the same particle during a calculation in 



which the Vh terms were included and equation (31) was used to calculate the smoothing 
lengths. A similar calculation was also performed, with the smoothing lengths being 
calculated according to the relation ([3^), and Njar = 6. The results from this calculation 
are presented in Fig. (14.c) 
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It is apparent from Fig. (14. b) that the use of equation (31) does lead to increased statis- 
tical fluctuations in the pressure forces experienced by the particles. These fluctuations, 
however, are generally not large relative to the absolute value of the overall pressure forces, 
and are found to decrease for increasing numbers of particles. The results shown in Fig. 
(14. c) show that it is possible to more or less remove these excess statistical fluctuations 
by simply increasing the number of particles contributing to the functional form used to 
calculate the smoothing lengths. The fact that the contribution of only a small number 
of particles {Nfar = 6) is required to reduce the pressure force fluctuations adds weight 
to the argument that the V/i terms should be regarded as essentially being correction 
terms. Up to now, however, we have not in general used equation (34) to calculate the 
smoothing lengths since the optimal method for implementing this relation has yet to 
be found. Finding the furthest Nfar nearest neighbours for each particle is a non trivial 
computational problem when efficiency is an important consideration. 



6.2 Numerical diffusion. 

Calculations of oscillating n = 3/2 polytropes were performed in order to examine whether 
the increased statistical fluctuations have the effect of increasing the numerical diffusivity 
of SPH. These used Ntol = 50 =b 1. At time t = 0, oscillations were initiated in an 
equilibrium polytrope, consisting of = 4945 particles, by imposing an homologous 
contraction. The results are presented in Fig. (15.) which, when moving from top to 
bottom in each panel, shows the time dependence of the thermal energy, kinetic energy, 
total energy, and gravitational potential energy. Fig. (15. a) shows the time sequence 
for a calculation in which the V/i terms were included and equation ( |3T|) was used to 
calculate the smoothing lengths, and Fig. (15. b) displays the oscillations for a similar 
calculation in which the V/i terms were neglected. In both of these calculations, a small 
amount (a = 0.1, /3 = 0.) of artificial viscosity was used in order to stabilise the time 
integration scheme, though the dissipated kinetic energy was not converted into thermal 
energy. Instead, the energy removed by the viscosity was simply allowed to leave the 
system, thus the polytropic equation of state was adopted throughout. 
In both calculations, the period of oscillation is accurate to within a few per cent of that 
expected from analytical considerations (see Cox 1980). Contrary to initial expectations, 
however, it is found that the oscillations damp out more quickly if the V/i terms are not 
included in the calculations, since the improved energy conservation gained by including 
the V/i terms more than off-sets any detrimental effects due to increased fiuctuations in 
the pressure forces. If the amplitude of the oscillations is measured by the height of the 
thermal energy curves in Figs. (15. a) and (15. b), then the oscillations damp by « 35% 
if the V/i terms are included, and by ~ 58% if the V/i terms are ignored, after thirteen 
oscillation periods. 



7 Discussion and conclusions. 

A reformulation of the SPH equations of motion that contain terms accounting for the 
local variability of the smoothing lengths has been presented. In a previous paper (Nelson 
& Papaloizou 1993) a set of conservative equations for a barotropic fiuid were derived. In 
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this current paper the derivation was generahsed to adiabatic fluids in which the entropy 
on each particle is conserved, and then extended in a self-consistent manner to incorporate 
the effects of viscous heating. Other non-adiabatic effects such as molecular cooling and 
radiative diffusion may be easily incorporated in a similar way, ensuring that the net rate 
of change of energy in the system due to such processes is equal to the time rate of change 
of the system Hamiltonian. 

The initial set of numerical tests presented here were designed to allow a direct compar- 
ison between a standard formulation of SPH employing spatially and temporally varying 
smoothing lengths, and our formulation of SPH in which the so-called V/i terms are in- 
cluded. It is apparent from the results of these calculations that the inclusion of the V/i 
terms has no detrimental effects on the ability of the algorithm to reproduce, with rea- 
sonable accuracy, the qualitative and quantitative features of known problems. In fact, 
owing to the improved energy conservation, it is found that certain features of the adia- 
batic collapse of a cold gas sphere are more accurately modeled when the V/i terms are 
included. 

The advantages of including the V/i terms in SPH are, however, most starkly brought to 
light by the numerical calculations specifically designed to test the conservation of energy. 
It is obvious from Figs. (12.) and (13.) that including the V/i terms leads to a dramatic 
improvement in conservation properties when the entropy equation ( p7| ) is integrated. 
Although total energy is formally conserved when the internal energy equation ( |33|) is 
used and the V/i terms are neglected in equations ( |3^ ) and (p3|), it is simple to show that 
entropy is no longer conserved when modeling adiabatic flows (e.g. Hernquist 1993). If the 
V/i terms are included, however, then equations (27) and (^) can be shown to be strictly 
equivalent, leading to the conservation of entropy and energy in both cases. Computational 
advantage may be obtained by choosing to integrate the entropy equation ( p7| ) since the 
use of velocities that are not properly time-centred in ( |27| ) still results in excellent energy 
conservation throughout the bulk of the flow. Test calculations summarised in Table 
3. imply that the problem of non-conservation is not due to gravitational, viscous, or 
discretisation effects, but is due to the time variability of the /ijS driving the evolution of 
the total energy. 

Up until now, discussions concerning V/i terms have tended to view the role of these 
additional terms as providing information about the local spatial variability of the smooth- 
ing lengths with respect to one another. However, once a functional form for the smoothing 
lengths has been defined in terms of inter-particle distances, as in equation (pO[), this in- 
terpretation appears not to be correct. Instead, it is apparent from equation (|3^ ) that the 
V/i terms provide information about the way VW{rij, hi) changes when we make changes 
to hi, independently of the way in which changes are made to the smoothing lengths of 
neighbouring particles. Thus, it is no longer obvious that one may justify neglecting V/i 
terms on the grounds that local spatial variations of the smoothing lengths are likely to 
be small because quantities do not vary greatly over the scale of the hiS. Instead, the 
V/i terms are accounting for the local time variability of the smoothing lengths, as the 
particles change their relative positions, and as such can only be neglected when the time 
scale over which the smoothing lengths vary appreciably is longer than both the dynamical 
time scale, and the time period over which the integration is performed. This is often not 
the case in calculations performed with SPH such as the collapse of gas clouds, or stellar 
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collisions. It is for the above reasons that the calculation performed by Hernquist (1993) in 
which each particle was given its own individual smoothing length, which was not allowed 
to vary in time, showed no substantial errors in energy conservation. V/i terms should 
not arise in such a situation since, in the Lagrangian sense, the smoothing lengths are not 
then functions of the particle positions/spatial coordinates. 

Finally, although it is relatively simple to derive a set of conservative equations from 
the Hamiltonian in SPH, taking into account the variability of the smoothing lengths, the 
practical problem of defining a functional form for the /ijS still remains. An important 
consideration is computational expense. If too many particles are allowed to contribute to 
the calculation of the smoothing lengths then the scheme becomes too costly in terms of 
CPU time. In order to avoid this problem, we have allowed only one particle to contribute 
to each hi, as in equation (31), for most of the calculations in this paper. At first glance 
this approach seems a little dangerous because additional fluctuations are then introduced 
into the pressure forces. In fact, by inspecting equation (|32| ) one might expect the V/i 
terms to be of the same order as the usual force term, in which case the fluctuations would 
be very large. However, the quantity dW/dh has a change of sign, with dW/dh < if 
< u/h < 1, and dW/dh > for 1 < u/h < 2. This leads to a process of self-cancellation, 
resulting in a much more diminished contribution from the V/i terms than one may have 
initially expected. A problem may still exist for particles that do not play the role of 
imax for other particles in the system, or for particles that play the role of imax for too 
many particles. In such a situation, the V/i contribution to the pressure force becomes 
either uni-directional, or else is of large magnitude due to the additive nature of these 
terms. For a purely random distribution of particles, it would be expected that such 
a situation might occur quite regularly. In practice, however, it seems that the system 
undergoes a process of relaxation in which correlations are set up between the particles, 
so that such occurrences are reduced to a few rogue incidents that have little effect on 
the overall evolution of the system. This is akin to the process described by Monaghan 
(1988) in order to explain why error estimates based on Monte-Carlo theory are larger 
that those observed in practice from SPH calculations. The test calculations presented in 
Section (6) show that in general these fluctuations are indeed small, and do not lead to 
an increase in the numerical diffusivity of SPH. In fact, the inclusion of V/i terms reduces 
the numerical diffusion in an oscillating polytrope due to improved energy conservation. 
However, we would still advise caution when using the functional form (|3l| ) for calculating 
the smoothing lengths, and suggest that any effects that may arise from its use be studied 
on a case by case basis. Should the increased pressure force fluctuations become a cause 
for concern, then one can easily increase the number of particles contributing to the /ijS 
through equation (p4[), though at the cost of increased CPU time and more complicated 
book-keeping. However, it seems from Fig. (14. c) that only a small number of particles 
(~ 6) are required to all but remove the excess fluctuations. 
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Tabic 1: Colliding Polytropcs. 

Run Vh terms Equation Ntol E- deviation (%) 

CI Included Entropy 50 0.8% 

C2 Included Entropy 50 ± 1 2.3% 

C3 Included Entropy 50 ± 2 3.2% 

C4 Excluded Entropy 50 ± 2 9.8% 

C5 Included Energy 50 ± 2 1.8% 

C6 Excluded Energy 50 ± 2 2.1% 
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Tabic 2: Expanding Polytropcs. 

Run Vh terms Equation Ntol E- deviation (%) 

El Included Entropy 50 0.01% 

E2 Included Entropy 50 ± 2 0.15% 

E^ Excluded Entropy 50 ± 2 1.48% 

EA Excluded Entropy 50 1.48% 

Ef> Included Energy 50 ± 2 0.13% 

EQ Excluded Energy 50 ± 2 0.13% 
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Table 3: Timing Tests for Colliding Polytropes: N=4170 



V/i terms /ij-eqn Mtol Time/t-step(s) 



Included 31 50 



95 



Included 31 50 ± 1 



81 



Included 31 50 ± 2 



78 



Excluded 31 50 ± 2 



72 



Included 34 50 ± 1 



95 



Included 34 50 ± 2 



87 
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Table Captions. 



Table 1. Results from colliding polytropes. Column(2) identifies whether the V/i terms 
were included, column(3) whether the entropy or internal energy equation was integrated. 
Column(4) gives the number of nearest neighbours within 2/i,, and column(5) the percent- 
age deviation from perfect energy conservation. 

Table 2. Results of timing tests performed for calculations of colliding polytropes. Col- 
umn (1) indicates whether V/i terms were included. Column (2) describes which equation, 
(|3l| ) or (^4[), was used to calculate smoothing lengths. The values taken for Mtol are given 
in column (3), and the CPU per time step, measured in seconds, is given in column (4). 

Table 3. Results from expanding polytropes. Column(2) identifies whether the Vh terms 
were included, column(3) whether the entropy or internal energy equation was integrated. 
Column(4) gives the number of nearest neighbours within 2/ij, Mtol^ and column(5) the 
percentage deviation from perfect energy conservation. 
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Figure Captions. 



Figure 1. (a) Density, (b) pressure, (c) velocity, (d) entropic profiles in a one-dimensional 
shock tube, using artificial viscosity defined by equation (p3|). Analytical solution is shown 
by dashed lines. In this case the V/i terms were neglected. 

Figure 2. (a) Density, (b) pressure, (c) velocity, (d) entropic profiles in a one-dimensional 
shock tube, using artificial viscosity defined by equation (^). In this case the V/i terms 
were included. 

Figure 3. Density profile during adiabatic collapse of a gas sphere, initially having an 
isothermal energy distribution and a 1/r density profile. Following Evrard (1988), density 
is measured in units of = SMt/^ttR"^ , where R is initial radius and Mt is total mass. 
Dimensionless time, normalised to free-fall time at outer radius, is shown in top right 
corner of each frame. In this case V/i terms were neglected. 

Figure 4. Thermal energy distribution, normalised to ti* = GMt/R, during adiabatic 
collapse of gaseous sphere. Dimensionless time is shown in top right corner of each frame. 
V/i terms were neglected. 

Figure 5. Pressure during adiabatic collapse of gas sphere, normalised to = p^u^. 
Dimensionless time is shown at top right corner of each frame. V/i terms were neglected. 

Figure 6. Radial velocity profiles during collapse of gas sphere described in text, nor- 
malised to = {GMt / R)^^^'^'^ ■ Dimensionless time is shown in top left corner of each 
frame. V/i terms were neglected. 

Figure 7. Same as Fig. 3, except that V/i terms were included. 
Figure 8. Same as Fig. 4, except that V/i terms were included. 
Figure 9. Same as Fig. 5, except that V/i terms were included. 
Figure 10. Same as Fig. 6, except that V/i terms were included. 

Figure 11. Head on collision between identical n = 3/2 polytropes, each containing 
2085 particles. Time is indicated at the top right corner of each panel, measured in 
computational units. Initial velocities and set-up is described in the text. 

Figure 12. Evolution of total energy in simulations of head on collisions between iden- 
tical n = 3/2 polytropes. The calculations represented in above figure correspond to the 
following runs listed in Table 1: CI - Solid line; C2 - Dashed line; C3 - Dot-dash-dotted 
line; C4 - Dotted line; C6 - Dash-dot-dot-dot-dashed line. Units are in computational 
units. 

Figure 13. Evolution of total energy in simulations of freely expanding polytropes in 
absence of gravity and viscosity. Units are in computational units. The calculations 
represented correspond to the following runs listed in Table 2: El - Solid line; E2 - 
Dashed line; E3 - Dash-dot-dashed line; E6 - Dotted line. 
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Figure 14. The graphs represent the pressure force, experienced as a function of time, 
by a single particle during free expansion of an n = 3/2 polytrope. (a) V/i terms are not 
included in equations of motion, (b) \7h terms are included and smoothing lengths are 



Figure 15. The graphs represent the various energies, as functions of time, during the 
normal mode oscillations of an n = 3/2 polytrope. In each panel, the line styles represent- 
ing each of the various energies are: dashed - thermal energy; dotted - kinetic energy; solid 
- total energy (including energy dissipated by viscosity); dash-dot-dashed - gravitational 
potential energy. Units are in computational units, (a). V/i terms included, (b). V/i 
terms neglected. 
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